{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "16ed848c-1758-4e62-bba6-20e82ef96dc2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.9/site-packages/numpy/ma/core.py:3380: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  _data[indx] = dval\n",
      "/opt/conda/lib/python3.9/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n",
      "/opt/conda/lib/python3.9/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n",
      "/opt/conda/lib/python3.9/site-packages/numpy/ma/core.py:3380: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  _data[indx] = dval\n",
      "/opt/conda/lib/python3.9/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAD8CAYAAABaU0PFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh60lEQVR4nO3df5Bc1XXg8e+xELZCSGSHAUsDDjhWJiHx2sQUDkVlaxzHFmidoLDJLmQrxomrZCdhK6QSYlinYieuXdhQ+eEUDlixKfDGwbgSkMlascza6SWuEuGXBELGAxgTrBGFCuPBjJkF/Tj7x7xxWq1uqef179ffT9XUdL++3X2Ouu/Tmfvuuy8yE0mSJEnL84pBByBJkiSNIgtpSZIkqQQLaUmSJKkEC2lJkiSpBAtpSZIkqQQLaUmSJKkEC2lJqpCIOC0i/ikiHomI3RHx28X210TEnRHxWPH71S2ef35EzETE4xFxZX+jl6TREq4jLUnVERFrgDWZ+UBEnAjcD2wE3gM8l5nXFAXyqzPzAw3PXQE8CrwD2APcC1ySmV/tYwqSNDIckZakCsnMpzPzgeL2C8AjwCRwIXBz0exmFovrRucAj2fmE5n5MvCZ4nmSpCaOG3QAzaxevTrf8IY3DDqMjn33u9/lhBNOGHQYXWEuw6cqeUC1crn//vufzcyJQccBEBGnA2cB/wKckplPw2KxHREnN3nKJPDNuvt7gLe2eO1NwCaAV73qVW953ete18XIB+fQoUO84hWjP8ZUlTzAXIZRVfIAePTRRzvaZw9lIX3KKadw3333DTqMjtVqNaanpwcdRleYy/CpSh5QrVwi4l8HHQNARHw/8PfA5Zn5nYho62lNtjWd/5eZm4HNAFNTUzkzM1M21KFSle9iVfIAcxlGVckDOt9nV+PPCUnS90TEShaL6E9n5m3F5meK+dNL86j3NXnqHuC0uvunAnt7GaskjTILaUmqkFgcev4k8Ehm/lndQ3cAlxa3LwU+1+Tp9wLrIuKMiDgeuLh4niSpCQtpSaqW84BfBX42InYWPxuAa4B3RMRjLK7KcQ1ARKyNiK0AmXkAuAzYxuJJip/NzN2DSEKSRsFQzpGWJJWTmV+h+VxngLc3ab8X2FB3fyuwtTfRSVK1OCItSZIklXDMQtqrZEmSJElHamdqxwHgd+uvkhURd7J4lawv1V0l60qg2VWyPkbdVbIi4o5hu0rWlh2zXLtthr1zC6xdvYor1k+x8azJQYclqQT7sySpX445Il31q2Rt2THLVbftYnZugQRm5xa46rZdbNkxO+jQJC2T/VmS1E/LOtmwX1fJmpiYoFarLSe00j5Se5GF/Ydfb2Bh/0E+8rkHWf38Yx299vz8fN/y6DVzGT5VyQO6l0sv+7MkSY3aLqT7fZWsfl0x57kvfL759v+XHV+1p0pX/jGX4VOVPKB7ufSyP0uS1KitVTuqfJWstatXLWu7pOFlf5Yk9VM7q3ZU+ipZV6yfYtXKFYdtW7VyBVesnxpQRJLKsj9LkvqpnRHpSl8la+NZk1x90Rs5fsXiP8Xk6lVcfdEbPctfGkH2Z0lSPx1zjvQ4XCVr41mT3HLPUwDc+r5zBxyNpE7YnyVJ/eKVDSVJkqQSLKQlSZKkEiykJUmSpBIspCVJkqQSLKQlSZKkEiykJUmSpBIspCVJkqQSLKQlSZKkEo55QRZJ0miJiBuBdwH7MvMni223AkvXSl8NzGXmm5s890ngBeAgcCAzz+5DyJI0kiykJal6bgKuAz61tCEz//PS7Yj4U+D5ozz/bZn5bM+ik6SKsJCWpIrJzLsi4vRmj0VEAP8J+Nm+BiVJFeQcaUkaLz8DPJOZj7V4PIEvRsT9EbGpj3FJ0shxRFqSxsslwC1Hefy8zNwbEScDd0bE1zLzrsZGRZG9CWBiYoJardaTYPttfn6+ErlUJQ8wl2FUlTy6wUJaksZERBwHXAS8pVWbzNxb/N4XEbcD5wBHFNKZuRnYDDA1NZXT09O9CLnvarUaVcilKnmAuQyjquTRDU7tkKTx8XPA1zJzT7MHI+KEiDhx6TbwTuDhPsYnSSPFQlqSKiYibgG2A1MRsSci3ls8dDEN0zoiYm1EbC3ungJ8JSIeBO4BPp+ZX+hX3JI0apzaIUkVk5mXtNj+nibb9gIbittPAG/qaXCSVCHHLKRd2F+SJEk6Ujsj0jfhwv6SJEnSYY5ZSLuwvyRJknSkTudIt7uwfwIfL5ZLamrQa5LOzS0AdPV9q7TOorkMn6rkAd3PpRf9WZKkRp0W0l1Z2B8Gvybp9TPbAZiePrdrr1mldRbNZfhUJQ/ofi696M+SJDUqvfxd3cL+t7ZqU7+wP7C0sL8kSZI08jpZR9qF/SVJkjS2jllIu7C/JEmSdKR2Vu1wYX9JkiSpgZcIlyRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJakiomIGyNiX0Q8XLftwxExGxE7i58NLZ57fkTMRMTjEXFl/6KWpNFjIS1J1XMTcH6T7X+emW8ufrY2PhgRK4CPARcAZwKXRMSZPY1UkkaYhbQkVUxm3gU8V+Kp5wCPZ+YTmfky8Bngwq4GJ0kVctygA5Ak9c1lEfFu4D7gdzPz2w2PTwLfrLu/B3hrsxeKiE3AJoCJiQlqtVr3ox2A+fn5SuRSlTzAXIZRVfLoBgtpSRoP1wMfAbL4/afArze0iSbPy2Yvlpmbgc0AU1NTOT093bVAB6lWq1GFXKqSB5jLMKpKHt3g1A5JGgOZ+UxmHszMQ8BfsziNo9Ee4LS6+6cCe/sRnySNomMW0p79LUmjLyLW1N39ReDhJs3uBdZFxBkRcTxwMXBHP+KTpFHUzoj0TXj2tySNjIi4BdgOTEXEnoh4L/AnEbErIh4C3gb8TtF2bURsBcjMA8BlwDbgEeCzmbl7IElI0gg45hzpzLwrIk4v8drfO/sbICKWzv7+aonXkiS1KTMvabL5ky3a7gU21N3fChwxOCJJOlInJxt27exvGPwZ4HNzCwBdfd8qndVqLsOnKnlA93PpRX+WJKlR2UK6q2d/w+DPAL9+ZjsA09Pndu01q3RWq7kMn6rkAd3PpRf9WZKkRqVW7fDsb0mSJI27UoW0Z39LkiRp3B1zakdx9vc0cFJE7AE+BExHxJtZnKrxJPC+ou1a4BOZuSEzD0TE0tnfK4AbPftbkiRJVdHOqh2e/S1JkiQ18MqGkiRJUgkW0pIkSVIJFtKSJElSCRbSkiRJUgkW0pIkSVIJFtKSJElSCRbSkiRJUgkW0pIkSVIJFtKSJElSCRbSkiRJUgkW0pIkSVIJFtKSJElSCRbSklQxEXFjROyLiIfrtl0bEV+LiIci4vaIWN3iuU9GxK6I2BkR9/UtaEkaQRbSklQ9NwHnN2y7E/jJzPx3wKPAVUd5/tsy882ZeXaP4pOkSrCQlqSKycy7gOcatn0xMw8Ud+8GTu17YJJUMccNOgBJUt/9OnBri8cS+GJEJPDxzNzcrFFEbAI2AUxMTFCr1XoRZ9/Nz89XIpeq5AHmMoyqkkc3HLOQjogbgXcB+zLzJ4tt1wI/D7wMfB34tcyca/LcJ4EXgIPAAQ8TStJgRcQHgQPAp1s0OS8z90bEycCdEfG1YoT7MEWBvRlgamoqp6enexVyX9VqNaqQS1XyAHMZRlXJoxvamdpxE861k6SRFxGXsjgw8l8yM5u1ycy9xe99wO3AOf2LUJJGyzELaefaSdLoi4jzgQ8Av5CZL7Zoc0JEnLh0G3gn8HCztpKk7syR7niuHQx+vt3c3AJAV9+3SnOIzGX4VCUP6H4uvejPoyQibgGmgZMiYg/wIRaPHL6SxekaAHdn5vsjYi3wiczcAJwC3F48fhzwt5n5hQGkIEkjoaNCultz7WDw8+2un9kOwPT0uV17zSrNITKX4VOVPKD7ufSiP4+SzLykyeZPtmi7F9hQ3H4CeFMPQ5OkSim9/J1z7SRJkjTOShXSzrWTJEnSuDtmIV3MtdsOTEXEnoh4L3AdcCKL0zV2RsQNRdu1EbG1eOopwFci4kHgHuDzzrWTJElSVRxzjrRz7SRJkqQjeYlwSZIkqQQLaUmSJKkEC2lJkiSpBAtpSZIkqQQLaUmSJKkEC2lJkiSpBAtpSZIkqQQLaUmSJKkEC2lJkiSpBAtpSZIkqQQLaUmSJKkEC2lJkiSpBAtpSZIkqQQLaUmSJKkEC2lJqpiIuDEi9kXEw3XbXhMRd0bEY8XvV7d47vkRMRMRj0fElf2LerC27JjlvGu+zHu+8F3Ou+bLbNkxO+iQpKFjPzmShbQkVc9NwPkN264EvpSZ64AvFfcPExErgI8BFwBnApdExJm9DXXwtuyY5arbdjE7twDA7NwCV922yyJBqmM/ae6YhbQjG5I0WjLzLuC5hs0XAjcXt28GNjZ56jnA45n5RGa+DHymeF6lXbtthoX9Bw/btrD/INdumxlQRNLwsZ80d1wbbW4CrgM+VbdtaWTjmqJAvhL4QP2T6kY23gHsAe6NiDsy86vdCHyYbdkxy7XbZpidW2Dy7i9zxfopNp41OeiwpKFjX+mrUzLzaYDMfDoiTm7SZhL4Zt39PcBbm71YRGwCNgFMTExQq9W6G20fLY2wNds+qnnNz8+PbOyNzGU4VLGfdMMxC+nMvCsiTm/YfCEwXdy+GajRUEhTN7IBEBFLIxuVLqSXDn0s/dW2dOgDsECQ6thXhlI02ZbNGmbmZmAzwNTUVE5PT/cwrN6avPvLTYuEydWrGNW8arXayMbeyFyGQxX7STeUnSN92MgG0O7IRuX/d/TQh9Qe+0rfPRMRawCK3/uatNkDnFZ3/1Rgbx9iG6gr1k+xauWKw7atWrmCK9ZPDSgiafjYT5prZ2pHWW2PbMDgDxPOFX9ldfq+VT30McqHoxpVJZdRz6OXfaVb/bli7gAuBa4pfn+uSZt7gXURcQYwC1wM/ErfIhyQpSMgv/93D/HywUNMrl7lNCOpgf2kubKF9DMRsaaYZ9eVkY1BHya8fmY7ANPT53b0OlU99DHKh6MaVSWXUc+jl32lW/15VEXELSxOvzspIvYAH2KxgP5sRLwXeAr45aLtWuATmbkhMw9ExGXANmAFcGNm7h5EDv228axJbrnnKebm5tj2gZ8ddDjSULKfHKns1I6lkQ1oY2QjIo5ncWTjjpLvNzI89CG1x77SO5l5SWauycyVmXlqZn4yM7+VmW/PzHXF7+eKtnszc0Pdc7dm5o9m5o9k5n8fXBaSNPzaWf7uFmA7MBURe4rRjGuAd0TEYyyuynFN0XZtRGwFyMwDwNLIxiPAZ8dhZGPjWZNcfdEbOX7F4j/t5OpVXH3RG8f+0IfUyL4iSRp17azacUmLh97epO1e4LCRDWBr6ehGlIc+pPbYVyRJo8wrG0qSJEklWEhLkiRJJVhIS5IkSSVYSEuSJEklWEhLkiRJJVhIS5IkSSVYSEuSJEklWEhLkiRJJVhIS5IkSSVYSEuSJEklWEhLkiRJJVhIS5IkSSVYSEuSJEklWEhLkiRJJVhIS5IkSSVYSEuSJEklWEhL0piIiKmI2Fn3852IuLyhzXREPF/X5g8HFK4kDb3jyj4xIqaAW+s2vR74w8z8i7o208DngG8Um27LzD8u+56SpPIycwZ4M0BErABmgdubNP3nzHxXH0OTpJFUupB2hyxJI+3twNcz818HHYgkjarShXQDd8iSNFouBm5p8di5EfEgsBf4vczc3dggIjYBmwAmJiao1Wq9irNv5uYWOHjwYCVymZ+fr0QeYC7Dpkr9pBu6VUh3tEOGwe+U5+YWALr2vlX7olWh8y+pSi5VyaMXfaXb/blqIuJ44BeAq5o8/ADww5k5HxEbgC3AusZGmbkZ2AwwNTWV09PTPYu3X66f2c7c3BxVyKVWq1UiDzCXYVOlftINHRfS3dghw+B3ytfPbAdgevrcrr1elb5oVej8S6qSS1Xy6EVf6XZ/rqALgAcy85nGBzLzO3W3t0bEX0XESZn5bF8jlKQR0I1VO466Q87M+eL2VmBlRJzUhfeUJJV3CS2OIkbEayMiitvnsPj/xLf6GJskjYxuTO046g4ZeCYz0x2yJA1eRHwf8A7gfXXb3g+QmTcAvwT8RkQcABaAizMzBxGrJA27jgppd8iSNFoy80Xghxq23VB3+zrgun7HJUmjqKNC2h2yJEmSxpVXNpQkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJKsJCWJEmSSrCQliRJkkqwkJYkSZJK6KiQjognI2JXROyMiPuaPB4R8ZcR8XhEPBQRP9XJ+0mSJEnD4rguvMbbMvPZFo9dAKwrft4KXF/8liQNQEQ8CbwAHAQOZObZDY8H8FFgA/Ai8J7MfKDfcUrSKOhGIX00FwKfyswE7o6I1RGxJjOf7vH7SpJacwBEkrqg00I6gS9GRAIfz8zNDY9PAt+su7+n2HZEIR0Rm4BNABMTE9RqtQ5DW565uQWArr3v3NwCBw8e7HsevTI/P28uQ6YqefSir3S7P48ZB0AkqU2dFtLnZebeiDgZuDMivpaZd9U9Hk2ek81eqCjCNwNMTU3l9PR0h6Etz/Uz2wGYnj63a683NzdHv/PolVqtZi5Dpip59KKvdLs/V0xXBkAGPfjRC1UaAKnKH9pgLsOmSv2kGzoqpDNzb/F7X0TcDpwD1BfSe4DT6u6fCuzt5D0lSR3pygDIoAc/eqFKAyBV+UMbzGXYVKmfdEPpVTsi4oSIOHHpNvBO4OGGZncA7y5W7/hp4HkPD0rS4NQPgABLAyD1HACRpDZ1svzdKcBXIuJB4B7g85n5hYh4f0S8v2izFXgCeBz4a+A3O4pWklSaAyCS1F2lp3Zk5hPAm5psv6HudgK/VfY9JElddQpw++IKdxwH/O3SAAh8b/+9lcWl7x5ncfm7XxtQrJI09Hq9/J0kaUg4ACJJ3eUlwiVJkqQSLKQlSZKkEiykJUmSpBIspCVJkqQSLKQlSZKkEiykJUmSpBIspCVJkqQSLKQlSZKkEiykJUmSpBIspCVJkqQSvES4KmnLjlmu3TbD3rkF1q5exRXrp9h41mTpdpIkSY0spFU5W3bMctVtu1jYfxCA2bkFrrptF8BhRXK77SRJkppxaocq59ptM98rjpcs7D/ItdtmSrWTJElqxhFpVc7euYW2trfbTpLUW07H06hyRFqVs3b1qra2t9tOktQ7S9PsZucWSP5tmt2WHbOl2kn9VLqQjojTIuKfIuKRiNgdEb/dpM10RDwfETuLnz/sLFzp2K5YP8WqlSsO27Zq5QquWD9Vqp0kqXecjqdR1snUjgPA72bmAxFxInB/RNyZmV9taPfPmfmuDt5HWpalw3y//3cP8fLBQ0y2OPzXbjtJUu84HU+jrHQhnZlPA08Xt1+IiEeASaCxkJb6buNZk9xyz1MA3Pq+cztuJ0nqjbWrVzHbpBhuNh2vnXZSP3XlZMOIOB04C/iXJg+fGxEPAnuB38vM3S1eYxOwCWBiYoJardaN0No2V3TObr3v3NwCBw8e7HsevTI/Pz9yubT6TBtz6fZn3y+j+Jk004u+Mqqfaa9FxGnAp4DXAoeAzZn50YY208DngG8Um27LzD/uY5gaM1esnzpsKVJoPR2vnXZSP3VcSEfE9wN/D1yemd9pePgB4Iczcz4iNgBbgHXNXiczNwObAaampnJ6errT0Jbl+pntAExPd2dU8vqZ7czNzdHvPHqlVquNXC6tPtPGXLr92ffLKH4mzfSir4zqZ9oHTsnT0HE6nkZZR4V0RKxksYj+dGbe1vh4fWGdmVsj4q8i4qTMfLaT95UkLZ9T8jSsnI6nUdXJqh0BfBJ4JDP/rEWb1xbtiIhzivf7Vtn3lCR1RztT8iLiHyPiJ/obmSSNjk5GpM8DfhXYFRE7i23/DXgdQGbeAPwS8BsRcQBYAC7OzOzgPSVJHerGlLxBn9fSC1U6t2UUz6Go+nktMJqfS6Mq9ZNu6GTVjq8AcYw21wHXlX0PSVJ3dWtK3qDPa+mFKp3bMornUFT9vBYYzc+lUZX6STd4ZUNJGhNOyZOk7urK8neSpJHglDxJ6iILafXVlh2zXLtthr1zC6wdk6WLxjFnDSen5ElSd1lIq2+27Jg9bDH92bkFrrptF0BlC8txzFmSpHHhHGn1zbXbZg67IhXAwv6DXLttZkAR9d445ixJ0rhwRFp9s7dYsqjd7VUwjjlLqpZxm542bvmqM45Iq2/Wrl61rO1VMI45S6qOpelps3MLJP82PW3LjtlBh9YT45avOmchrb65Yv0Uq1auOGzbqpUruGL91IAi6r1xzFlSdYzb9LRxy1edG/upHUuHcGbnFjh+xSvYsmPWQzg9svTv+vt/9xAvHzzE5BgcMhvHnIfBsy+8xDe/vcAZV37eQ7NSB8Ztetq45avOjXUh3biiwssHD3H5rTv5o3/YzYd+/if8j7cHNp41yS33PAXAre8bvatSlTGOOQ9C/R/F9VwpRSpv7epVR/Sppe1VNG75qnNjPbWj2SEcgG+/uN85UdIIqZ/X2IyHZqVyxm162rjlq86NdSF9tEM1/scrjY5WfxTX89CstHwbz5rk6oveyPErFsuFydWruPqiN1b26M645avOjfXUjlaHcJb4H680Gtrpqx6alcoZt+lp45avOjPWI9Jv+7GJo14r9wdXrVz2a27ZMct513yZf/nGc8x8+xCnX/l5zrvmy04TkRp0s6+s/r6j99Vgsb9LktRNY1tI/8GWXXz67qfIo7SZW9jPH2zZ1fZrtpqn6TqU0uG62Vf+YMsuvv3i/qO2SeDTdz+1rP4sSdKxjNXUji07ZvnwHbuZWzj6f7r1/ubup/ibu5/i1d+38qgreWzZMcvvfHYn2aIyX9h/kA/fsdt5VhLwR/+wu+Wc5oX9B/mdz+4EWq+yUaYvJ+33Z0mS2jGUhfST31k8zDtMvv3ifi6/dSeX37qz9GvMLezvWV4WBuq2MsVqt2TScX87mm7051aOf+0b3tL1F5UkDaWOCumIOB/4KLAC+ERmXtPweBSPbwBeBN6TmQ908p5qrpeFwfd8obt/BPTrj6Wm79Mkl37E0/X36PJnIkmS2ld6jnRErAA+BlwAnAlcEhFnNjS7AFhX/GwCri/7fpIkSdIw6WRE+hzg8cx8AiAiPgNcCHy1rs2FwKcyM4G7I2J1RKzJzKeP9sKnvrCP//HPf9VBaJLG2coVr+DAoUMtz1nopXf3/y0lSQPSyaodk8A36+7vKbYttw0AEbEpIu6LiPs6iKm01a8Kfuw1K1j9qqMtiCdp2K1+VfAjPxhMvdr+LEnqrU5GpJv9D9U4/tNOm8WNmZuBzQCvXLMuP/Azv9lBaO1bOklvfXGS3pvozklWJxy/gtd9/yEeeW4AQ2LSiPnx1wRPvhAs7D9U+jUa+zJ0rz8vyxOX9+d9JEkDF1ny2GdEnAt8ODPXF/evAsjMq+vafByoZeYtxf0ZYPpYUzteuWZdrrn0L0rF1Y7lrHCx3P+E61+7Vqsx94PrBrbygTTslvrL6ucfY3p6uqP+1o5+FNVP33w5Lz392FAOhffqBPFe77MlqVc63Wd3UkgfBzwKvB2YBe4FfiUzd9e1+Q/AZSzulN8K/GVmnnOs156amsqZmZlScQ2TWq3G9PR0T99jkEuUaTz0Y2nFfvSVfomI+zPz7EHH0ag4QfxR4B0sTrO7F7gkM79a12YD8F/5t332RzPzrcd6bQtpSaOq00K69NSOzDwQEZcB21gc3bgxM3dHxPuLx28AtrK4Q36cxdGNXyv7fmpu41mTfVk7ukqFTlVyqUoe6htPEJekBp2eIN7ROtKZuZXFYrl+2w11txP4reW+7qOPPjpfTAMZdScBzw46iC4xl+FTlTygWrlMDTqAFpqd/N042tzqBPEjCumI2MTisqa8YtUP8O4nHu1qsJLUDwee39fR84fyyobAzDAeGl2uiLivCnmAuQyjquQB1ctl0DG00LMTxCPivpdefL4yn18VvotVyQPMZRhVJQ/ofJ/dyfJ3kqTRsQc4re7+qcDeEm0kSQULaUkaD/cC6yLijIg4HrgYuKOhzR3Au2PRTwPPH2t+tCSNs2Gd2rF50AF0SVXyAHMZRlXJA8yl53p8gvhQ5lxSVXKpSh5gLsOoKnlAh7mUXv5OkiRJGmdO7ZAkSZJKsJCWJEmSShhoIR0RvxwRuyPiUEScXbf99IhYiIidxc8NdY+9JSJ2RcTjEfGXxSVtB65VLsVjVxXxzkTE+rrtQ5lLvYj4cETM1n0WG+oea5rXsIqI84tYH4+IKwcdz3JFxJPF92Xn0nI9EfGaiLgzIh4rfr960HE2ExE3RsS+iHi4blvL2If5u9Uil8r0k6Nxnz2cudSr0nfRffbguM9eRi6ZObAf4MdZvHhBDTi7bvvpwMMtnnMPcC6L653+I3DBIHNoI5czgQeBVwJnAF8HVgxzLg15fRj4vSbbW+Y1jD8snlz1deD1wPFF7GcOOq5l5vAkcFLDtj8BrixuXwn8z0HH2SL2fw/8VH2/bhX7sH+3WuRSiX7SRu7us4cwl4a8KvFddJ898NjdZ7eZy0BHpDPzkcxs+wqGEbEG+IHM3J6LGX8K2Nir+JbjKLlcCHwmM1/KzG+weDb8OcOcS5ua5jXgmI7me5dHzsyXgaXLI4+6C4Gbi9s3M6Tfocy8C3iuYXOr2If6u9Uil1aGOpflcp89nLm0adS+i+6zB8h9dvu5DPMc6TMiYkdE/N+I+Jli2ySLFwxYsnT52mHW6pK7o5TLZRHxUHF4ZOlQTqu8htWoxdtMAl+MiPtj8fLMAKdksc5v8fvkgUW3fK1iH9XPqgr9pBPus4dHFb6LoxZvM+6zh1tX+knP15GOiP8DvLbJQx/MzM+1eNrTwOsy81sR8RZgS0T8BMu4fG0vlMylVcwDzaXe0fICrgc+wmJsHwH+FPh1hij+No1avM2cl5l7I+Jk4M6I+NqgA+qRUfysqtJP3Ge7zx4WoxZvM+6zh1fX+knPC+nM/LkSz3kJeKm4fX9EfB34URb/Mji1rmlfL19bJhdaX3J3oLnUazeviPhr4H8Xd0ftUsKjFu8RMnNv8XtfRNzO4uGmZyJiTWY+XRx63jfQIJenVewj91ll5jNLt0e8n7jPdp89LEYt3iO4zx5e3dxnD+XUjoiYiIgVxe3XA+uAJ4pDCS9ExE8XZ0u/G2g1qjAs7gAujohXRsQZLOZyz6jkUnSWJb8ILJ312jSvfse3DO1cHnloRcQJEXHi0m3gnSx+FncAlxbNLmUIv0NH0Sr2UftuVamflOI+e3hU6LvoPnv4uM9uZsBnUv4ii9X/S8AzwLZi+38EdrN45uQDwM/XPefsIuGvA9dRXJ1x0D+tcike+2AR7wx1Z3kPay4Nef0vYBfwUPEFW3OsvIb1h8VLHz9axPzBQcezzNhfX/SHB4u+8cFi+w8BXwIeK36/ZtCxtoj/FhYP/+8v+sl7jxb7MH+3WuRSmX5yjNzdZw9hLg15Vea76D57oPG7z24zFy8RLkmSJJUwlFM7JEmSpGFnIS1JkiSVYCEtSZIklWAhLUmSJJVgIS1JkiSVYCEtSZIklWAhLUmSJJXw/wFU0025SA7u/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "\n",
    "# 定义一个函数，用于生成正弦波\n",
    "def getSin(amp, freq, phase, sampleList):\n",
    "    return amp * np.sin(-2 * math.pi * freq * sampleList + phase)\n",
    "\n",
    "\n",
    "# 定义一个函数，用于生成余弦波\n",
    "def getCos(amp, freq, phase, sampleList):\n",
    "    return amp * np.cos(-2 * math.pi * freq * sampleList + phase)\n",
    "\n",
    "\n",
    "# 定义去噪函数，将振幅小于阈值的噪声去除\n",
    "def denoise(arr, thresh):\n",
    "    mask = arr > thresh\n",
    "    return mask * arr\n",
    "\n",
    "\n",
    "# 设置采样率为 3000\n",
    "srate = 3000\n",
    "# 在 0 到 1 的范围内生成等间距的采样点\n",
    "t = np.linspace(0, 1, srate)\n",
    "\n",
    "# 生成四个正弦和余弦波\n",
    "signals = []\n",
    "signals.append(getSin(amp=1.5, freq=30, phase=0, sampleList=t))\n",
    "signals.append(getCos(amp=3, freq=5, phase=0, sampleList=t))\n",
    "signals.append(getSin(amp=10, freq=100, phase=0, sampleList=t))\n",
    "signals.append(getCos(amp=20, freq=120, phase=0, sampleList=t))\n",
    "\n",
    "# 将四个正弦和余弦波相加，得到混合信号\n",
    "mixed_signal = np.zeros_like(signals[0])\n",
    "for signal in signals:\n",
    "    mixed_signal += signal\n",
    "\n",
    "# 使用 numpy 的 fft 函数计算傅里叶系数\n",
    "fCoefs = np.fft.fft(mixed_signal, srate)\n",
    "\n",
    "# 计算振幅列表\n",
    "amp_list = np.zeros_like(fCoefs)\n",
    "for i in range(len(fCoefs)):\n",
    "    amp_list[i] = 2 * np.abs(fCoefs[i] / srate)\n",
    "\n",
    "# 将振幅数组和频率数组进行平移，使 0 Hz 位于中心\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1 / srate)\n",
    "amp_shifted = np.fft.fftshift(amp_list)\n",
    "freq_shift = np.fft.fftshift(freqs)\n",
    "\n",
    "# 创建两个子图，分别用于显示原始振幅频谱和处理后的振幅频谱\n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "# 绘制原始振幅频谱\n",
    "ax[0].stem(freq_shift, amp_shifted)\n",
    "ax[0].set_xlim([-150, 150])\n",
    "ax[0].set_ylim([-0.5, 21])\n",
    "ax[0].grid()\n",
    "\n",
    "# 对振幅频谱进行去噪处理\n",
    "for i in range(len(freq_shift)):\n",
    "    if freq_shift[i] > 110 or freq_shift[i] < -110:\n",
    "        amp_shifted[i] = 0\n",
    "amp_shifted = denoise(amp_shifted, 1)\n",
    "\n",
    "\n",
    "# 绘制处理后的振幅频谱\n",
    "ax[1].stem(freq_shift, amp_shifted)\n",
    "# 设置第二个子图的坐标轴范围和刻度\n",
    "ax[1].set_xlim([-150, 150])\n",
    "ax[1].set_ylim([0, 20])\n",
    "ax[1].set_yticks(np.arange(0, 21, 2.5))\n",
    "ax[1].set_yticklabels(np.arange(0, 21, 2.5))\n",
    "ax[1].grid()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# 定义一个函数，用于生成正弦波\n",
    "def getSin(amp, freq, phase, sampleList):\n",
    "    return amp * np.sin(-2 * math.pi * freq * sampleList + phase)\n",
    "\n",
    "# 定义一个函数，用于生成余弦波"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f3d0b931-ea99-4db9-88d9-67cb35dcbf8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAD8CAYAAABaU0PFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkXklEQVR4nO3df5SddX3g8feHIWD4oVEJSAIWXOlssVapHCyH056xisGsLamru9A9VavbqC17Sk9LhbVHbT1badn+sAcLRmXBrUVcC5GuqZHVzqLnBEFIIEQcQKSYCUsKYYAhA0kmn/3j3sGbm3szd577+5n365w5c+/zfJ97P5/c+33yme/zfZ4nMhNJkiRJC3NYvwOQJEmShpGFtCRJklSAhbQkSZJUgIW0JEmSVICFtCRJklSAhbQkSZJUgIW0JJVIRJwcEf8cEfdFxLaI+N3q8pdFxC0R8UD190ubbH9eRExExIMRcWlvo5ek4RJeR1qSyiMiTgROzMy7IuJY4E5gDfBeYFdmXl4tkF+amR+u23YEuB84F9gO3AFcmJnf72EKkjQ0HJGWpBLJzEcz867q42eA+4CVwPnAddVm11EpruudBTyYmQ9l5h7gS9XtJEkNHN7vABpZtmxZvvrVr+53GG179tlnOfroo/sdRkeYy+ApSx5QrlzuvPPOxzNzeb/jAIiIU4AzgO8CJ2Tmo1AptiPi+AabrAR+XPN8O/DGJq+9FlgL8KIXvegNr3zlKzsYef/s37+fww4b/jGmsuQB5jKIypIHwP3339/WPnsgC+kTTjiB733ve/0Oo23j4+OMjY31O4yOMJfBU5Y8oFy5RMS/9DsGgIg4BvgH4OLMfDoiWtqswbKG8/8ycx2wDmB0dDQnJiaKhjpQyvJdLEseYC6DqCx5QPv77HL8OSFJekFELKFSRH8xM2+sLn6sOn96bh71zgabbgdOrnl+ErCjm7FK0jCzkJakEonK0PPngfsy8y9rVt0MvKf6+D3AVxtsfgdwWkScGhFHABdUt5MkNWAhLUnlcg7wG8AvR8SW6s9q4HLg3Ih4gMpVOS4HiIgVEbEBIDP3ARcBG6mcpPjlzNzWjyQkaRgM5BxpSVIxmfkdGs91Bnhzg/Y7gNU1zzcAG7oTnSSViyPSkiRJUgHzFtLeJUuSJEk6WCtTO/YBv197l6yIuIXKXbK+WXOXrEuBRnfJ+jQ1d8mKiJsH7S5Z6zdPcsXGCXZMzbBi2VIuWTXKmjNW9jssSQXYnyVJvTLviHTZ75K1fvMkl924lcmpGRKYnJrhshu3sn7zZL9Dk7RA9mdJUi8t6GTDXt0la/ny5YyPjy8ktMI+Mb6bmb0H3m9gZu8sn/jq3Sx76oG2Xnt6erpneXSbuQyesuQBnculm/1ZkqR6LRfSvb5LVq/umLPr619rvPy5bPuuPWW684+5DJ6y5AGdy6Wb/VmSpHotXbWjzHfJWrFs6YKWSxpc9mdJUi+1ctWOUt8l65JVoyxdMnLAsqVLRrhk1WifIpJUlP1ZktRLrYxIl/ouWWvOWMkn3/Fajhip/FOsXLaUT77jtZ7lLw0h+7MkqZfmnSO9GO6SteaMlVx/+yMA3PCBs/scjaR22J8lSb3inQ0lSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpALmvSGLJGm4RMQ1wNuBnZn5s9VlNwBz90pfBkxl5usbbPsw8AwwC+zLzDN7ELIkDSULaUkqn2uBK4EvzC3IzP849zgi/gJ46hDbvykzH+9adJJUEhbSklQymXlrRJzSaF1EBPAfgF/uaVCSVELOkZakxeUXgccy84Em6xP4RkTcGRFrexiXJA0dR6QlaXG5ELj+EOvPycwdEXE8cEtE/CAzb61vVC2y1wIsX76c8fHxrgTba9PT06XIpSx5gLkMorLk0QkW0pK0SETE4cA7gDc0a5OZO6q/d0bETcBZwEGFdGauA9YBjI6O5tjYWDdC7rnx8XHKkEtZ8gBzGURlyaMTnNohSYvHW4AfZOb2Risj4uiIOHbuMfBW4N4exidJQ8VCWpJKJiKuBzYBoxGxPSLeX111AXXTOiJiRURsqD49AfhORNwN3A58LTO/3qu4JWnYOLVDkkomMy9ssvy9DZbtAFZXHz8EvK6rwUlSicxbSHthf0mSJOlgrYxIX4sX9pckSZIOMG8h7YX9JUmSpIO1O0e61Qv7J/CZ6uWSGur3NUmnpmYAOvq+ZbrOorkMnrLkAZ3PpRv9WZKkeu0W0h25sD/0/5qkV01sAmBs7OyOvWaZrrNoLoOnLHlA53PpRn+WJKle4cvf1VzY/4ZmbWov7A/MXdhfkiRJGnrtXEfaC/tLkiRp0Zq3kPbC/pIkSdLBWrlqhxf2lyRJkup4i3BJkiSpAAtpSZIkqQALaUmSJKkAC2lJkiSpAAtpSZIkqQALaUmSJKkAC2lJkiSpAAtpSZIkqQALaUmSJKkAC2lJkiSpAAtpSSqZiLgmInZGxL01yz4eEZMRsaX6s7rJtudFxEREPBgRl/YuakkaPhbSklQ+1wLnNVj+V5n5+urPhvqVETECfBp4G3A6cGFEnN7VSCVpiFlIS1LJZOatwK4Cm54FPJiZD2XmHuBLwPkdDU6SSuTwfgcgSeqZiyLi3cD3gN/PzCfr1q8EflzzfDvwxkYvFBFrgbUAy5cvZ3x8vPPR9sH09HQpcilLHmAug6gseXSChbQkLQ5XAZ8Asvr7L4D31bWJBttloxfLzHXAOoDR0dEcGxvrWKD9ND4+ThlyKUseYC6DqCx5dIJTOyRpEcjMxzJzNjP3A5+lMo2j3nbg5JrnJwE7ehGfJA2jeQtpz/6WpOEXESfWPP014N4Gze4ATouIUyPiCOAC4OZexCdJw6iVEelr8exvSRoaEXE9sAkYjYjtEfF+4M8jYmtE3AO8Cfi9atsVEbEBIDP3ARcBG4H7gC9n5ra+JCFJQ2DeOdKZeWtEnFLgtV84+xsgIubO/v5+gdeSJLUoMy9ssPjzTdruAFbXPN8AHDQ4Ikk6WDsnG3bs7G/o/xngU1MzAB193zKd1Woug6cseUDnc+lGf5YkqV7RQrqjZ39D/88Av2piEwBjY2d37DXLdFaruQyesuQBnc+lG/1ZkqR6ha7a4dnfkiRJWuwKFdKe/S1JkqTFbt6pHdWzv8eA4yJiO/AxYCwiXk9lqsbDwAeqbVcAn8vM1Zm5LyLmzv4eAa7x7G9JkiSVRStX7fDsb0mSJKmOdzaUJEmSCrCQliRJkgqwkJYkSZIKsJCWJEmSCrCQliRJkgqwkJYkSZIKsJCWJEmSCrCQliRJkgqwkJYkSZIKsJCWJEmSCrCQliRJkgqwkJYkSZIKsJCWpJKJiGsiYmdE3Fuz7IqI+EFE3BMRN0XEsibbPhwRWyNiS0R8r2dBS9IQspCWpPK5FjivbtktwM9m5s8B9wOXHWL7N2Xm6zPzzC7FJ0mlYCEtSSWTmbcCu+qWfSMz91Wf3gac1PPAJKlkDu93AJKknnsfcEOTdQl8IyIS+ExmrmvUKCLWAmsBli9fzvj4eDfi7Lnp6elS5FKWPMBcBlFZ8uiEeQvpiLgGeDuwMzN/trrsCuBXgD3AD4HfzMypBts+DDwDzAL7PEwoSf0VER8B9gFfbNLknMzcERHHA7dExA+qI9wHqBbY6wBGR0dzbGysWyH31Pj4OGXIpSx5gLkMorLk0QmtTO24FufaSdLQi4j3UBkY+U+ZmY3aZOaO6u+dwE3AWb2LUJKGy7yFtHPtJGn4RcR5wIeBX83M3U3aHB0Rx849Bt4K3NuorSSpM3Ok255rB/2fbzc1NQPQ0fct0xwicxk8ZckDOp9LN/rzMImI64Ex4LiI2A58jMqRwyOpTNcAuC0zPxgRK4DPZeZq4ATgpur6w4G/z8yv9yEFSRoKbRXSnZprB/2fb3fVxCYAxsbO7thrlmkOkbkMnrLkAZ3PpRv9eZhk5oUNFn++SdsdwOrq44eA13UxNEkqlcKXv3OunSRJkhazQoW0c+0kSZK02M1bSFfn2m0CRiNie0S8H7gSOJbKdI0tEXF1te2KiNhQ3fQE4DsRcTdwO/A159pJkiSpLOadI+1cO0mSJOlg3iJckiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSSiYiromInRFxb82yl0XELRHxQPX3S5tse15ETETEgxFxae+i7q/1myc55/Jv8d6vP8s5l3+L9Zsn+x2SNHDsJwezkJak8rkWOK9u2aXANzPzNOCb1ecHiIgR4NPA24DTgQsj4vTuhtp/6zdPctmNW5mcmgFgcmqGy27capEg1bCfNDZvIe3IhiQNl8y8FdhVt/h84Lrq4+uANQ02PQt4MDMfysw9wJeq25XaFRsnmNk7e8Cymb2zXLFxok8RSYPHftLY4S20uRa4EvhCzbK5kY3LqwXypcCHazeqGdk4F9gO3BERN2fm9zsR+CBbv3mSKzZOMDk1w8rbvsUlq0ZZc8bKfoclDRz7Sk+dkJmPAmTmoxFxfIM2K4Ef1zzfDryx0YtFxFpgLcDy5csZHx/vbLQ9NDfC1mj5sOY1PT09tLHXM5fBUMZ+0gnzFtKZeWtEnFK3+HxgrPr4OmCcukKampENgIiYG9kodSE9d+hj7q+2uUMfgAWCVMO+MpCiwbJs1DAz1wHrAEZHR3NsbKyLYXXXytu+1bBIWLlsKcOa1/j4+NDGXs9cBkMZ+0knFJ0jfcDIBtDqyEbp/3f00IfUGvtKzz0WEScCVH/vbNBmO3ByzfOTgB09iK2vLlk1ytIlIwcsW7pkhEtWjfYpImnw2E8aa2VqR1Etj2xA/w8TTlX/ymr3fct66GOYD0fVK0suw55HN/tKp/pzydwMvAe4vPr7qw3a3AGcFhGnApPABcCv9yzCPpk7AvKHX7mHPbP7WblsqdOMpDr2k8aKFtKPRcSJ1Xl2HRnZ6PdhwqsmNgEwNnZ2W69T1kMfw3w4ql5Zchn2PLrZVzrVn4dVRFxPZfrdcRGxHfgYlQL6yxHxfuAR4F3VtiuAz2Xm6szcFxEXARuBEeCazNzWjxx6bc0ZK7n+9keYmppi44d/ud/hSAPJfnKwolM75kY2oIWRjYg4gsrIxs0F329oeOhDao19pXsy88LMPDEzl2TmSZn5+cx8IjPfnJmnVX/vqrbdkZmra7bdkJk/nZn/JjP/W/+ykKTB18rl764HNgGjEbG9OppxOXBuRDxA5aocl1fbroiIDQCZuQ+YG9m4D/jyYhjZWHPGSj75jtdyxEjln3blsqV88h2vXfSHPqR69hVJ0rBr5aodFzZZ9eYGbXcAB4xsABsKRzekPPQhtca+IkkaZt7ZUJIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSJEkqwEJakiRJKsBCWpIkSSrAQlqSFomIGI2ILTU/T0fExXVtxiLiqZo2H+1TuJI08A4vumFEjAI31Cx6FfDRzPzrmjZjwFeBH1UX3ZiZf1L0PSVJxWXmBPB6gIgYASaBmxo0/XZmvr2HoUnSUCpcSLtDlqSh9mbgh5n5L/0ORJKGVeFCuo47ZEkaLhcA1zdZd3ZE3A3sAP4gM7fVN4iItcBagOXLlzM+Pt6tOHtmamqG2dnZUuQyPT1dijzAXAZNmfpJJ3SqkG5rhwz93ylPTc0AdOx9y/ZFK0Pnn1OWXMqSRzf6Sqf7c9lExBHArwKXNVh9F/BTmTkdEauB9cBp9Y0ycx2wDmB0dDTHxsa6Fm+vXDWxiampKcqQy/j4eCnyAHMZNGXqJ53QdiHdiR0y9H+nfNXEJgDGxs7u2OuV6YtWhs4/pyy5lCWPbvSVTvfnEnobcFdmPla/IjOfrnm8ISL+NiKOy8zHexqhJA2BTly145A75Mycrj7eACyJiOM68J6SpOIupMlRxIh4RURE9fFZVP6feKKHsUnS0OjE1I5D7pCBxzIz3SFLUv9FxFHAucAHapZ9ECAzrwbeCXwoIvYBM8AFmZn9iFWSBl1bhbQ7ZEkaLpm5G3h53bKrax5fCVzZ67gkaRi1VUi7Q5YkSdJi5Z0NJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpAIspCVJkqQCLKQlSZKkAiykJUmSpALaKqQj4uGI2BoRWyLiew3WR0T8TUQ8GBH3RMTPt/N+kiRJ0qA4vAOv8abMfLzJurcBp1V/3ghcVf0tSeqDiHgYeAaYBfZl5pl16wP4FLAa2A28NzPv6nWckjQMOlFIH8r5wBcyM4HbImJZRJyYmY92+X0lSc05ACJJHdBuIZ3ANyIigc9k5rq69SuBH9c8315ddlAhHRFrgbUAy5cvZ3x8vM3QFmZqagagY+87NTXD7Oxsz/PolunpaXMZMGXJoxt9pdP9eZFxAESSWtRuIX1OZu6IiOOBWyLiB5l5a836aLBNNnqhahG+DmB0dDTHxsbaDG1hrprYBMDY2Nkde72pqSl6nUe3jI+Pm8uAKUse3egrne7PJdORAZB+D350Q5kGQMryhzaYy6ApUz/phLYK6czcUf29MyJuAs4Cagvp7cDJNc9PAna0856SpLZ0ZACk34Mf3VCmAZCy/KEN5jJoytRPOqHwVTsi4uiIOHbuMfBW4N66ZjcD765eveMXgKc8PChJ/VM7AALMDYDUcgBEklrUzuXvTgC+ExF3A7cDX8vMr0fEByPig9U2G4CHgAeBzwK/3Va0kqTCHACRpM4qPLUjMx8CXtdg+dU1jxP4naLvIUnqqBOAmypXuONw4O/nBkDghf33BiqXvnuQyuXvfrNPsUrSwOv25e8kSQPCARBJ6ixvES5JkiQVYCEtSZIkFWAhLUmSJBVgIS1JkiQVYCEtSZIkFWAhLUmSJBVgIS1JkiQVYCEtSZIkFWAhLUmSJBVgIS1JkiQV4C3CVUrrN09yxcYJdkzNsGLZUi5ZNcqaM1YWbidJklTPQlqls37zJJfduJWZvbMATE7NcNmNWwEOKJJbbSdJktSIUztUOldsnHihOJ4zs3eWKzZOFGonSZLUiCPSKp0dUzMtLW+1nSSpu5yOp2HliLRKZ8WypS0tb7WdJKl75qbZTU7NkPxkmt36zZOF2km9VLiQjoiTI+KfI+K+iNgWEb/boM1YRDwVEVuqPx9tL1xpfpesGmXpkpEDli1dMsIlq0YLtZMkdY/T8TTM2pnasQ/4/cy8KyKOBe6MiFsy8/t17b6dmW9v432kBZk7zPeHX7mHPbP7Wdnk8F+r7SRJ3eN0PA2zwoV0Zj4KPFp9/ExE3AesBOoLaann1pyxkutvfwSAGz5wdtvtJEndsWLZUiYbFMONpuO10k7qpY6cbBgRpwBnAN9tsPrsiLgb2AH8QWZua/Iaa4G1AMuXL2d8fLwTobVsqto5O/W+U1MzzM7O9jyPbpmenh66XJp9pvW5dPqz75Vh/Ewa6UZfGdbPtNsi4mTgC8ArgP3Ausz8VF2bMeCrwI+qi27MzD/pYZhaZC5ZNXrApUih+XS8VtpJvdR2IR0RxwD/AFycmU/Xrb4L+KnMnI6I1cB64LRGr5OZ64B1AKOjozk2NtZuaAty1cQmAMbGOjMqedXEJqampuh1Ht0yPj4+dLk0+0zrc+n0Z98rw/iZNNKNvjKsn2kPOCVPA8fpeBpmbRXSEbGEShH9xcy8sX59bWGdmRsi4m8j4rjMfLyd95UkLZxT8jSonI6nYdXOVTsC+DxwX2b+ZZM2r6i2IyLOqr7fE0XfU5LUGa1MyYuIf4qI1/Q2MkkaHu2MSJ8D/AawNSK2VJf9V+CVAJl5NfBO4EMRsQ+YAS7IzGzjPSVJberElLx+n9fSDWU6t2UYz6Eo+3ktMJyfS70y9ZNOaOeqHd8BYp42VwJXFn0PSVJndWpKXr/Pa+mGMp3bMoznUJT9vBYYzs+lXpn6SSd4Z0NJWiSckidJndWRy99JkoaCU/IkqYMspNVT6zdPcsXGCXZMzbBikVy6aDHmrMHklDxJ6iwLafXM+s2TB1xMf3Jqhstu3ApQ2sJyMeYsSdJi4Rxp9cwVGycOuCMVwMzeWa7YONGniLpvMeYsSdJi4Yi0emZH9ZJFrS4vg8WYs6RyWWzT0xZbvmqPI9LqmRXLli5oeRksxpwllcfc9LTJqRmSn0xPW795st+hdcViy1fts5BWz1yyapSlS0YOWLZ0yQiXrBrtU0TdtxhzllQei2162mLLV+1b9FM75g7hTE7NcMTIYazfPOkhnC6Z+3f9w6/cw57Z/axcBIfMFmPOg+DxZ57nx0/OcOqlX/PQrNSGxTY9bbHlq/Yt6kK6/ooKe2b3c/ENW/jjf9zGx37lNf7H2wVrzljJ9bc/AsANHxi+u1IVsRhz7ofaP4preaUUqbgVy5Ye1KfmlpfRYstX7VvUUzsaHcIBeHL3XudESUOkdl5jIx6alYpZbNPTFlu+at+iLqQPdajG/3il4dHsj+JaHpqVFm7NGSv55DteyxEjlXJh5bKlfPIdry3t0Z3Flq/at6indjQ7hDPH/3il4dBKX/XQrFTMYpuettjyVXsW9Yj0m/7t8kPeK/clS5cs+DXXb57knMu/xXd/tIuJJ/dzyqVf45zLv+U0EalOJ/vKsqMO3VeDSn+XJKmTFm0h/Ufrt/LF2x4hD9FmamYvf7R+a8uv2WyeptehlA7Uyb7yR+u38uTuvYdsk8AXb3tkQf1ZkqT5LKqpHes3T/Lxm7cxNXPo/3Rr/d1tj/B3tz3CS49acsgreazfPMnvfXkL2aQyn9k7y8dv3uY8Kwn443/c1nRO88zeWX7vy1uA5lfZKNKXk9b7syRJrRjIQvrhpyuHeQfJk7v3cvENW7j4hi2FX2NqZm/X8rIwUKcVKVY7JZO2+9uhdKI/N3PEK179ho6/qCRpILVVSEfEecCngBHgc5l5ed36qK5fDewG3puZd7Xznmqsm4XBC77e2T8CevXHUsP3aZBLL+Lp+Ht0+DORJEmtKzxHOiJGgE8DbwNOBy6MiNPrmr0NOK36sxa4quj7SZIkSYOknRHps4AHM/MhgIj4EnA+8P2aNucDX8jMBG6LiGURcWJmPnqoFz7pmZ386bf/to3QJC1mS0YOY9/+/U3PWeimd/f+LSVJfdJOIb0S+HHN8+3AG1tosxI4qJCOiLVURq151dEvbiOsYpa9KHjFUYfx/3bvZ+q5PvzvK6kjKn05gBH7c7dEHANsA44C2BQBL385PPkkPfnrJQKWLoWZmQPf7+ijYffuyrIIOOooOPLIA7e7+mp45zu7H6O0mE1Pw+mnV/rjnMze7SPmU7N/2AWvI+JxKuekf4jMryzkpdoppBtdgrn+X6eVNpWFmeuAdQBHnnhafvgXf7uN0Fo3d5LequpJeq+jMydZHX3ECK88Zj/37RqAL4w04H7mZcHDzwQze/cXfo36vgyd688L8tDFvXmffsqcJuK3gC8DsQTgmWcq6+JQV+fvoOeeg+OOg8cfrzwfGYE9ew6MYWamUuCPVG/5/PrXW0RLvXDMMfC5z8G73vWTwvm55yq/e7WPmE91//Dsrl37XwpPA1sWWkQDRBb8yyAizgY+npmrqs8vA8jMT9a0+QwwnpnXV59PAGPzTe048sTT8sT3/HWhuFqxkCtcLPQ/4drXHh8fZ+olp/XtygfSoJvrL8ueeoCxsbG2+lsrelFUP3rdxTz/6AMD8j/FgTp+gnjEtcAvPRojp56Q+2kyTtI1s1E5zeewTHYd9WJevvvpg2KYjcN4/vAj2HfYCOf+56v412Ne1tMYpcXsv3/trzjrx/dy5L49LH92il7vI+YzG4fxXCbHkk8Cr2Ge+rSRdgrpw4H7gTcDk8AdwK9n5raaNv8OuIjKTvmNwN9k5lnzvfbo6GhOTEwUimuQjI+PMzY21tX36OclyrQ49OLSir3oK70SEXdm5pn9jqNe9QTx+4FzqUyzuwO4MDO/X9NmNfBf+Mk++1OZWT9lr/ZFXwrctx9O6GLohzRz+JHMHnYYx+xpfpv4mSUv4uNv+QD/6+fO7WFkkl783DTf/OwHePnup/odSlPPRnBs5vvI/B9Fti88tSMz90XERcBGKqMb12Tmtoj4YHX91cAGKjvkB6mMbvxm0fdTY2vOWNmTa0eXqdApSy5lyUM90/kTxDOfJOLehBP6NQR/5P697Dpq2SEL6SNn93D2E9/nzG8/0MPIJAHMLhnIW5a8YKQynnxj0e3byi4zN1AplmuXXV3zOIHfWejr3n///dPVaSDD7jjg8X4H0SHmMnjKkgeUK5fRfgfQRMdPEP8QXPoJOOlo+nN3rwSeSjjm6V3saRJDAjs4jCe2fpeLXnJ8bwOUFrkLZ55h7Omppv2z3xLYCZwCfwZ8sMhrDGJeABODeGh0oSLie2XIA8xlEJUlDyhfLv2OoYnOniAO64E/BWaehZEXVbbt9cD0zHG5/0XV990NLG0Qw8xP7d83+1szT+dvzTz9djJvPdQLluW7WJY8wFwGUUt5RBxP5YjXEir98hh6v4+Yz8zx5JHArxPx9/PtHxoZ1EJaktRZ24GTa56fBOwo0GbOlVTOdOdJOPJoeAZ4aWdCbckeKv8pz13qZQnwRF0Me+ZirPocEa8l8/nehCgtai/sI6qep7f7iPnsAZ7eDS8/qvKHeKH9Q+E7G0qShsodwGkRcWpEHAFcANxc1+Zm4N1R8QvAUw3nR0esAd4CLAOWvbgyKDNCZfQ6evADMA0cUbNshMqIV9a0+VdgpuZnCfDRBfybSSqibh9R/enlPqKVfci/AjPPV/4YL7x/GNQR6XX9DqBDypIHmMsgKkseYC5d19ETxDPXU5naAcBLItZW7wVQBuYxeMxl8Bw6j7p9xCA7qc39V+HL30mSJEmLmVM7JEmSpAIspCVJkqQC+lpIR8S7ImJbROyPiDNrlp8SETMRsaX6c3XNujdExNaIeDAi/qZ6S9u+a5ZLdd1l1XgnImJVzfKBzKVWRHw8IiZrPovVNesa5jWoIuK8aqwPRsSl/Y5noSLi4er3ZcvcJdYi4mURcUtEPFD9PUhnRL8gIq6JiJ0RcW/NsqaxD/J3q0kupeknh+I+ezBzqVWm76L77P5xn72AXDKzbz/Az1C5ecE4cGbN8lOAe5tscztwNpWzLv8JeFs/c2ghl9OBu4EjgVOBHwIjg5xLXV4fB/6gwfKmeQ3iD5WTq34IvIrKmf53A6f3O64F5vAwcFzdsj8HLq0+vhT4s37H2ST2XwJ+vrZfN4t90L9bTXIpRT9pIXf32QOYS11epfguus/ue+zus1vMpa8j0pl5X2a2fAfDiDgReHFmbspKxl8A1nQrvoU4RC7nA1/KzOcz80dUzoY/a5BzaVHDvPoc06G8cHvkzNwDzN0eedidD1xXfXwdA/odyspF7nfVLW4W+0B/t5rk0sxA57JQ7rMHM5cWDdt30X12H7nPbj2XQZ4jfWpEbI6I/xsRv1hdtpLKDQPmzN2+dpA1u+XuMOVyUUTcUz08Mncop1leg2rY4m0kgW9ExJ0Rsba67ISsXue3+nuY7oHcLPZh/azK0E/a4T57cJThuzhs8TbiPnuwdaSfdP060hHxf4BXNFj1kcz8apPNHgVemZlPRMQbgPUR8RoWcPvabiiYS7OY+5pLrUPlBVwFfIJKbJ8A/gJ4HwMUf4uGLd5GzsnMHVG57eotEfGDfgfUJcP4WZWln7jPdp89KIYt3kbcZw+ujvWTrhfSmfmWAts8T+VWkmTmnRHxQ+CnqfxlcFJN00PdvrbjiuRC81vu9jWXWq3mFRGfBf539elCbiU8CIYt3oNk5o7q750RcROVw02PRcSJmflo9dDzzr4GuTDNYh+6zyozH5t7POT9xH22++xBMWzxHsR99uDq5D57IKd2RMTyiBipPn4VcBrwUPVQwjMR8QvVs6XfDTQbVRgUNwMXRMSREXEqlVxuH5Zcqp1lzq8Bc2e9Nsyr1/EtQCu3Rx5YEXF0RBw79xh4K5XP4mbgPdVm72EAv0OH0Cz2YftulamfFOI+e3CU6LvoPnvwuM9upM9nUv4aler/eeAxYGN1+b8HtlE5c/Iu4FdqtjmzmvAPgSup3p2x3z/Ncqmu+0g13glqzvIe1Fzq8vqfwFbgnuoX7MT58hrUHyq3Pr6/GvNH+h3PAmN/VbU/3F3tGx+pLn858E3ggervl/U71ibxX0/l8P/eaj95/6FiH+TvVpNcStNP5sndffYA5lKXV2m+i+6z+xq/++wWc/EW4ZIkSVIBAzm1Q5IkSRp0FtKSJElSARbSkiRJUgEW0pIkSVIBFtKSJElSARbSkiRJUgEW0pIkSVIB/x8T9iYckXlvdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def getCos(amp, freq, phase, sampleList):\n",
    "    return amp * np.cos(-2 * math.pi * freq * sampleList + phase)\n",
    "\n",
    "# 定义去噪函数，将振幅小于阈值的噪声去除\n",
    "def denoise(arr, thresh):\n",
    "    mask = arr > thresh\n",
    "    return mask * arr\n",
    "\n",
    "# 设置采样率为 3000\n",
    "srate = 3000\n",
    "# 在 0 到 1 的范围内生成等间距的采样点\n",
    "t = np.linspace(0, 1, srate)\n",
    "\n",
    "# 生成四个正弦和余弦波\n",
    "signals = []\n",
    "signals.append(getSin(amp=1.5, freq=30, phase=0, sampleList=t))\n",
    "signals.append(getCos(amp=3, freq=5, phase=0, sampleList=t))\n",
    "signals.append(getSin(amp=10, freq=100, phase=0, sampleList=t))\n",
    "signals.append(getCos(amp=20, freq=120, phase=0, sampleList=t))\n",
    "\n",
    "# 将四个正弦和余弦波相加，得到混合信号\n",
    "mixed_signal = np.zeros_like(signals[0])\n",
    "for signal in signals:\n",
    "    mixed_signal += signal\n",
    "\n",
    "# 使用 numpy 的 fft 函数计算傅里叶系数\n",
    "fCoefs = np.fft.fft(mixed_signal)\n",
    "\n",
    "# 计算振幅列表\n",
    "amp_list = np.abs(fCoefs) / len(mixed_signal) * 2\n",
    "\n",
    "# 将振幅数组和频率数组进行平移，使 0 Hz 位于中心\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1 / srate)\n",
    "amp_shifted = np.fft.fftshift(amp_list)\n",
    "freq_shift = np.fft.fftshift(freqs)\n",
    "\n",
    "# 创建两个子图，分别用于显示原始振幅频谱和处理后的振幅频谱\n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "# 绘制原始振幅频谱\n",
    "ax[0].stem(freq_shift, amp_shifted)\n",
    "ax[0].set_xlim([-150, 150])\n",
    "ax[0].set_ylim([-0.5, 21])\n",
    "ax[0].grid()\n",
    "\n",
    "# 对振幅频谱进行去噪处理\n",
    "for i in range(len(freq_shift)):\n",
    "    if freq_shift[i] > 110 or freq_shift[i] < -110:\n",
    "        amp_shifted[i] = 0\n",
    "amp_shifted = denoise(amp_shifted, 1)\n",
    "\n",
    "# 绘制处理后的振幅频谱\n",
    "ax[1].stem(freq_shift, amp_shifted)\n",
    "# 设置第二个子图的坐标轴范围\n",
    "ax[1].set_xlim([-150, 150])\n",
    "ax[1].set_ylim([0, 20])\n",
    "ax[1].set_yticks(np.arange(0, 21, 2.5))\n",
    "ax[1].grid()\n",
    "#在第二个子图中标记去除的频率点\n",
    "for i in range(len(freq_shift)):\n",
    "      if freq_shift[i] > 110 or freq_shift[i] < -110:\n",
    "         ax[1].annotate('X', xy=(freq_shift[i], amp_shifted[i]), ha='center', va='center', color='r', fontsize=12)\n",
    "\n",
    "#显示图像\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "144f38a3-0c61-447f-9c03-1975553e8bc2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
